/**
  CNOK project, Anyang Normal University, IMP-CAS
  \class TAEikonalPhase
  \brief compute eikonal phase \chi=\chi_N + \chi_C, S=exp(i*\chi), where \chi_N
  represents nuclear phase, and \chi_C Coulomb phase.
  This is a static class.
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2020/07/25
  \date Last modified: 2020/09/10 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#include <cstring>
#include "config.h"
#ifdef CATCH_CNOK
#include "catch.h"
#endif
#include <fstream>
#include <unistd.h>
#include <libgen.h> // basename, dirname
#include <iostream>
#include "TAAME.h"
#include "TAFun.h"
#include "TAEikonalPhase.h"
#include "TAFNN.h"
#include "TABessel.h"
#include "TAIntegrate.h"
#include "TAMath.h"
#include "TAFourier.h"
#include "TAInterpolate.h"
#include "TAAME.h"
#include "TAParaManager.h"
#include "TAYaml.h"

using std::cout;
using std::endl;
using std::ifstream;
using std::ofstream;

#define polint TAInterpolateD::PolyInter
#define bessj0 TABessel::BesselJ0
#define romb TAIntegrate<>::Romberg
#define simp TAIntegrate<>::Simpson
#define endwith TAMath::endwith

// fine structure constant ~1/137, e^2/hbarc
static const double ALPHA = TAMath::FineStructureConstant();
static const double HBARC = TAMath::hbarc();
static const double u = TAMath::uMeV(); // atomic mass unit in MeV
static const double FOURPI = 4.*TAMath::Pi();
static const cdouble I(0., 1.); // imaginary number unit I
static const TAAME *AME = TAAME::Instance();
// static const TAParaManager *PARA = TAParaManager::Instance();

// the maximum momentum transfer in fm^-1, q=p/hbar*sin(\theta)
// e.g. for E=240 MeV/nucleon, QMAX=6 corresponds to scattering angle > 25 degree.
const double TAEikonalPhase::QMAX = 6.;

// alphaP, alphaT: proj(targ) nucleon size factor: rho_n(r)=exp(-r^2/alphap^2)
TAEikonalPhase::TAEikonalPhase(int zP, int aP, int zT, int aT, double ek,
  double alphaP, double alphaT) : fZP(zP), fAP(aP), fZT(zT), fAT(aT), fEk(ek),
    fAlphaP(alphaP), fAlphaT(alphaT){
  fPM = AME->massNUC(fZP, fAP);
  fTM = AME->massNUC(fZT, fAT);
  fMu = fPM*fTM/(fPM+fTM);
  if(zP < 0 || aP < 0 || zT < 0 || aT < 0 || ek < 0. || aP < zP || aT < zT)
    TAException::Error("TAEikonalPhase", "TAEikonalPhase: abnormal input.");

  // whether to use pauli-blocking or not //
  fIsPauli = false;
  const YAML::Node &pauli = (*PARA->GetPara())["isPauli"];
  if(pauli) fIsPauli = pauli.as<bool>();
  int nP = aP - zP, nT = aT - zT;
  // fNN(q)=kNN/4pi*sigmaNN(i+alphaNN)*exp(-betaNN*q^2)
  fBeta = TAMath::EkPerUToBeta(ek);
  fAlphaNN = TAFNN::GetAlphaNN(ek, zP, nP, zT, nT);
  fBetaNN  = TAFNN::GetBetaNN(ek, zP, nP, zT, nT);
  fSigmaNN = TAFNN::GetSigmaNN(ek, zP, nP, zT, nT, fIsPauli);
  cout << "Nucleon-nucleon scattering data: " << endl;
  cout << "Ek: " << ek << " MeV/nucleon" << endl;
  cout << "Isospin averaged NN-cross section sigmaNN: " << fSigmaNN*10. << " mb" << endl;
  cout << "Imaginary coefficient of the f_NN amplitude: " << fAlphaNN << endl;
  cout << "Exponential coefficient of the f_NN amplitude: " << fBetaNN << endl;

  fNRP = fNRT = 0;
  fRP = fRhoP = fRT = fRhoT = nullptr;
  fQ = fRhoQ = nullptr;

  fIsOpt = false; fOptN = 0;
  fOptR = fV = fW = nullptr;
} // end the constructor

TAEikonalPhase::~TAEikonalPhase(){
  if(fRP) delete [] fRP;
  if(fRhoP) delete [] fRhoP;
  if(fRT) delete [] fRT;
  if(fRhoT) delete [] fRhoT;
  if(fQ) delete [] fQ;
  if(fRhoQ) delete [] fRhoQ;

  if(fOptR) delete [] fOptR;
  if(fV) delete [] fV;
  if(fW) delete [] fW;
} // end the destructor

void TAEikonalPhase::SetProjectileDensity(const string &file){
  ReadDensityFile(file, fNRP, fRP, fRhoP, fAP);
} // end member function SetProjectileDensity
void TAEikonalPhase::SetTargetDensity(const string &file){
  ReadDensityFile(file, fNRT, fRT, fRhoT, fAT);
} // end member function SetTargetDensity
void TAEikonalPhase::ReadDensityFile(const string &file, int &n, double *&r, double *&rho, int a){
  // density is a delta function for empty file, so its 2D Fourier transform is a
  // constant, A=1, stored in n, and returned over by TAFourier::Fourier2D(...)
  if(r || rho) TAException::Error("TAEikonalPhase", "ReadDensityFile: r or rho is not nullptr");
  if("" == file){ n = a; return; } // single-nucleon, whose Fourier transform is a constant, i.e. A

  const char *f = file.c_str();
  if(access(f, F_OK))
    TAException::Error("TAEikonalPhase", "ReadDensityFile: %s does not exist.", f);
  ifstream ff(f);
  // get Nof lines in file //
  string line; int nl = 0; n = 0;
  while(getline(ff, line)) nl++;
  ff.clear(); ff.seekg(0);
  r = new double[nl]{}, rho = new double[nl]{};
  if(endwith(file, ".den")){
    getline(ff, line); // skip comments
    while(ff >> r[n] >> rho[n]) n++;
  } // end if
  else if(endwith(file, ".rho")){
    double rn, rp, rx;
    for(int i = 3; i--;) getline(ff, line); // skip comments
    while(ff >> r[n] >> rp >> rn >> rho[n] >> rx) n++;
  } // end else if
  else TAException::Error("TAEikonalPhase", "ReadDensityFile: %s: unknown file suffix.", f);
  if(!ff.eof()) TAException::Error("TAEikonalPhase", "ReadDensityFile: %s: format error", f);
  // normalize to nucleon number
  const double cc = TAFourier::Fourier2D(0., n, r, rho); // \int{rho(r)d^3r}
  for(int i = n; i--;) rho[i] *= a/cc;
  // cout << "a=" << a << " cc=" << cc << endl; getchar(); // DEBUG
} // end member function ReadDensityFile

// set optical potentials
// file format: r V W: V_opt(r)=V(r)+iW(r)
void TAEikonalPhase::SetOpticalPotential(const string &opt){
  fIsOpt = true;
  const char *f = opt.c_str();
  if(access(f, F_OK))
    TAException::Error("TAEikonalPhase", "SetOpticalPotential: %s does not exist.", f);
  ifstream ff(f);
  // get Nof lines in file //
  string line; int nl = 0;
  while(getline(ff, line)) nl++;
  ff.clear(); ff.seekg(0);
  fOptR = new double[nl]{}, fV = new double[nl]{}, fW = new double[nl]{};
  getline(ff, line); // skip comments
  while(ff >> fOptR[fOptN] >> fV[fOptN] >> fW[fOptN]) fOptN++;
  if(!ff.eof()) TAException::Error("TAEikonalPhase", "SetOpticalPotential: %s: format error", f);
} // end member function SetOpticalPotential

/// \retval eikonal nucclear phase = 1/KNN*\int_0^\infty{q*rhoP*rhoT*fNN*J0(qb)}
/// \param b: impact parameter
double TAEikonalPhase::GetPhaseNI(double b){
  // using optical potential //
  if(fIsOpt){
    if(0 == fOptN) TAException::Error("TAEikonalPhase", "GetPhaseNI: %s: fOptN = 0");
    const auto &vw = [=](double z){
      return polint(fOptR, fW, fOptN, 6, sqrt(b*b+z*z));
    }; // end lambda expression
    return -2./(HBARC*fBeta)*romb(vw, 0., 30.);
  } // end if

  if(!fNRP || !fNRT) TAException::Error("TAEikonalPhase",
    "GetPhaseNI: Densities may not be assigned yet.");
  // compute integral \int_0^\infty{q*rhoP*rhoT*fNN*J0(qb)}
  static const int NQ = 10000; // Nof grid of Q mesh
  // calculate the density's 2D Fourier transform once and for all //
  if(!fRhoQ){
    fQ = new double[NQ];
    fRhoQ = new double[NQ];
    const double dq = QMAX / (NQ-1);
    char name1[128], name2[128], name3[128], nc[8] = "v"; // DEBUG
    if(1 != fAP) sprintf(nc, "c"); // DEBUG
    sprintf(name1, "vis/%s.txt", nc); // DEBUG
    sprintf(name2, "vis/t.txt"); // DEBUG
    sprintf(name3, "vis/%st.txt", nc); // DEBUG
    cout << endl << "Calculating 2D Fourier transform of nuclear densities rho(q)..." << endl;
    cout << "rho_P(q) is stored in \033[32;1m" << name1 << "\033[0m." << endl;
    cout << "rho_T(q) is stored in \033[32;1m" << name2 << "\033[0m." << endl;
    cout << "q*rhoqP*rhoqT*exp(-fBetaNN*q2) is ";
    cout << "stored in \033[32;1m" << name3 << "\033[0m." << endl;
    ofstream fout1(name1), fout2(name2), fout3(name3); // DEBUG
    for(int i = 0; i < NQ; i++){
      const double q = fQ[i] = dq*i;
      const double q2 = q*q;
      // note that rho(q) could be minus for large q //
      const double rhoqP = TAFourier::Fourier2D(q, fNRP, fRP, fRhoP)*exp(-sqr(q*fAlphaP)/4.);
      const double rhoqT = TAFourier::Fourier2D(q, fNRT, fRT, fRhoT)*exp(-sqr(q*fAlphaT)/4.);
      fRhoQ[i] = q*rhoqP*rhoqT*exp(-fBetaNN*q2);
      fout1 << q << " " << rhoqP << endl; // DEBUG
      fout2 << q << " " << rhoqT << endl; // DEBUG
      fout3 << q << " " << fRhoQ[i] << endl; // DEBUG
      // cout << q << " QDESNP=" << rhoqP << endl; // DEBUG
      // cout << q << " QDESNT=" << rhoqT << endl; // DEBUG
      // cout << q << " fRHOQ=" << fRhoQ[i] << endl; // DEBUG
      // cout << q << " EF=" << exp(-fBetaNN*q2) << endl; // DEBUG
      // getchar(); // DEBUG
    } // end for over i
    cout << "...DONE. rho(q)'s can be inspected via ";
    cout << "`cd vis/` and `root -l rhoq.C`" << endl;
    cout << "(provided ROOT is installed)." << endl << endl;
    fout1.close(); fout2.close(); fout3.close();
  } // end if
  // calculate the integral 1/KNN*\int_0^\infty{q*rhoP*rhoT*fNN*J0(qb)} //
  double pshift[NQ]{};
  for(int i = 0; i < NQ; i++){
    // b=0.1; // DEBUG
    pshift[i] = fRhoQ[i]*bessj0(fQ[i]*b);
    // cout << "b=" << b << " q=" << fQ[i] << endl; // DEBUG
    // cout << "fRHOQ=" << fRhoQ[i] << endl; // DEBUG
    // cout << "bessj0(fQ[i]*b)=" << bessj0(fQ[i]*b) << endl; // DEBUG
    // getchar(); // DEBUG
  } // end for over i
  // const auto &ph = [=](double q){
  //   static const int NP = 4; // Nof points in fRhoQ interpolation
  //   return polint(fQ, fRhoQ, NQ, NP, q)*TABessel::BesselJ0(q*b);
  // }; // end of the lambda expression
  // cout << "DUM=" << simp(NQ, fQ, pshift) << endl; getchar(); // DEBUG
  return simp(NQ, fQ, pshift)*fSigmaNN/FOURPI;
  // return romb(ph, 0., QMAX, 1.e-3)*fSigmaNN/FOURPI;
} // end member function GetPhaseNI

double TAEikonalPhase::GetPhaseNR(double b){
  // using optical potential //
  if(fIsOpt){
    if(0 == fOptN) TAException::Error("TAEikonalPhase", "GetPhaseNR: %s: fOptN = 0");
    const auto &vv = [=](double z){
      return polint(fOptR, fV, fOptN, 6, sqrt(b*b+z*z));
    }; // end lambda expression
    return -2./(HBARC*fBeta)*romb(vv, 0., 30.);
  } // end if

  return GetPhaseNI(b)*fAlphaNN;
} // end member function GetPhaseNR

cdouble TAEikonalPhase::GetPhaseN(double b){
  return GetPhaseNI(b)*(I+fAlphaNN);
} // end member function GetPhaseN

/// \retval eikonal Coulumb phase = 2\eta*ln(kb)
/// \param b: impact parameter
double TAEikonalPhase::GetPhaseC(double b) const{
  if(0. == b) b = 1.e-4;
  // Sommerfeld parameter eta = zp*zt*e^2/(hbar*v)=zp*zt/beta * ALPHA
  const double eta = fZP*fZT/fBeta * ALPHA; // ALPHA: fine structure constant
  const double k = fMu*fBeta/HBARC; // wave number of AP-AT scattering system, in fm^(-1)
  // cout << "fZP: " << fZP << " fZT: " << fZT << endl; // DEBUG
  // cout << "eta: " << eta << " beta: " << fBeta << " ALPHA: " << ALPHA << endl; // DEBUG
  // cout << "k: " << k << " fMu: " << fMu << " HBARC: " << HBARC << endl; // DEBUG
  return 2.*eta*log(k*b); // 2.*\eta*ln(k*b)
  // return 0.;
} // end member function GetPhaseC

/// \retval \chi_C+\chi_N
/// \param b: impact parameter
cdouble TAEikonalPhase::GetPhase(double b){
  return GetPhaseN(b) + GetPhaseC(b);
} // end memeber function GetPhase

#ifdef CATCH_CNOK
TEST_CASE("template class TAIntegrate", "[quad]"){
  SECTION("polynomial tests"){
    class QuadTest : public TAFun<double>{
    public:
      double operator()(double x) const override{
        static const int n = 7;
        return (n+1)*pow(x, n);
      } // end member function
    }; // end of the class
    QuadTest quad;
    double t = TAIntegrate<>::Trapezoid(quad, 0., 1.);
    double s = TAIntegrate<>::Simpson(quad, 0., 1.);
    double r = TAIntegrate<>::Romberg(quad, 0., 1.);
    CHECK(t == Approx(1.).epsilon(1.e-4));
    CHECK(s == Approx(1.).epsilon(1.e-6));
    CHECK(r == Approx(1.).epsilon(1.e-8));
  } // end section0
  SECTION("trigonometric function tests"){
    class QuadTest : public TAFun<double>{
    public:
      double operator()(double x) const override{
        return sin(x);
      } // end member function
    }; // end of the class
    QuadTest quad;
    const double b = TAMath::Pi()/2.;
    double t = TAIntegrate<>::Trapezoid(quad, 0., b);
    double s = TAIntegrate<>::Simpson(quad, 0., b);
    double r = TAIntegrate<>::Romberg(quad, 0., b);
    CHECK(t == Approx(1.).epsilon(1.e-4));
    CHECK(s == Approx(1.).epsilon(1.e-6));
    CHECK(r == Approx(1.).epsilon(1.e-6));
  } // end section1
} // end TEST_CASE
#endif
